Myxozoan survey of thicklip grey mullet Chelon labrosus reinforces successful radiation of Myxobolus in mugiliform hosts

A myxozoan survey was performed on specimens of thicklip grey mullet Chelon labrosus (Risso) captured from the Douro River estuary, northern Portugal. Eleven new species, all belonging to the genus Myxobolus Bütschli, 1882 (M. abdominalis n. sp., M. aestuarium n. sp., M. caudalis n. sp., M. chelonari n. sp., M. cucurbitiformis n. sp., M. douroensis n. sp., M. intestinicola n. sp., M. invictus n. sp., M. labicola n. sp., M. peritonaei n. sp., and M. pinnula n. sp.) are described based on microscopic and molecular data, confirming the known high radiation of these myxozoans in mullets. Additionally, Myxobolus pupkoi Gupta et al., 2022 is reported for the first time from C. labrosus, bringing forth a novel case of morphological plasticity between geographic isolates. We consider that molecular-based comparisons are imperative for the description of mugiliform-infecting Myxobolus, with distance estimation further matching two of the novel Myxobolus spp. with sphaeractinomyxon types previously reported from another Portuguese estuary. This finding supports sphaeractinomyxon as specific life cycle counterparts of Myxobolus that infect mullets. Phylogenetic analyses of 18S rDNA retrieved a monophyletic clade of mugiliform-infecting myxobolids comprising well-supported lineages of species parasitizing mullets from the genera Chelon, Mugil, Crenimugil, and Planiliza. The existence of more than one Chelon- and Planiliza-infecting lineage reveals that myxobolids parasitized members of these genera multiple times during their evolution. Lastly, the elevated number of unmatched sphaeractinomyxon sequences included in the Chelon-infecting lineages clearly shows that Myxobolus diversity hosted by this genus remains underrated.


Introduction
Mugilids (Teleostei, Mugilidae) are highly opportunistic and diverse fishes that inhabit a wide range of habitats worldwide; because they have a high tolerance to varying abiotic conditions like salinity, temperature, sedimentary regimes, turbidity, and dissolved oxygen [46], mullets can be found in coastal areas, estuaries, lakes, and lagoons, in both tropical and temperate regions. This ubiquitous nature and extreme adaptability grant them excellent survival skills and allows them to perform important roles in ecosystems; they are key players in the flow of matter and energy from lower to higher trophic levels [1,5,25]. It also leaves them vulnerable and exposed to parasitic infections and disease.
Myxozoa Grassé, 1970 is a subphylum of cnidarians that evolved to become obligate endoparasites of invertebrates and vertebrates [36]. Myxozoan parasites are also ubiquitous/widespread in aquatic environments, existing both in fresh and saltwater habitats [30,36]. Myxozoan species are traditionally described using morphological characters as the main taxonomic criteria. However, descriptions based solely on morphology have been proved to be unreliable due to the extreme diversity of this group and high plasticity of myxospores between and within genera [2,3,17,18]. In fact, the molecular analyses of several gene markers have revealed discrepancies between morphology and molecular-based classifications, showing that morphological characters are not enough and may even be misleading when it comes to genus and specieslevel classification [2,19,28,39,40]. Therefore, a more integrative approach is required for performing reliable species identifications, which should be based not only on morphology, but also DNA sequencing, host identification, habitat, tissue tropism, histopathology, and other details of the sporogonic stages and their development [2,19].
The present study aimed to acknowledge the myxozoan diversity infecting C. labrosus in another northern Portuguese estuary (Douro River), contributing to the general knowledge of myxozoan morphological and molecular data, as well as life cycle clarification.

Materials and methods
Fish sampling and myxozoan survey Thirteen specimens of thicklip grey mullet C. labrosus (Risso, 1927) were acquired fresh from commercial stocks captured by local fishermen in the Douro River estuary, near the city of Porto (41°08 0 N8°39 0 W), during March 2022. Examination of the fish's external (fins, eyes, skin, and scales) and internal organs and tissues (brain, gills, heart, gonads, liver, stomach, intestines, spleen, gallbladder, swim bladder, urinary bladder, kidney, and muscle) was performed both macroscopically and microscopically, using an Olympus BX41 microscope (Olympus, Coimbra, Portugal) for the detection of myxospores, plasmodia, cysts, or other developmental stages of Myxozoa. Cysts and infected tissues were examined and photographed using an Olympus DP72 digital camera (Olympus). Measurements were taken from fresh mature myxospores, following the guidelines of Lom and Arthur [29]. All measurements are given in micrometers (lm), and include mean ± standard deviation, range of variation (minimum and maximum), and number of myxospores measured (n). The prevalence (%) of infection was determined based on matching of the 18S rDNA sequences obtained for each case isolate in study.

DNA extraction, amplification, and sequencing
After being photographed, cysts and infected tissue samples were separately preserved in absolute ethanol at 4°C. DNA was extracted with the extraction kit "GenElute TM Mammalian Genomic DNA Miniprep" (Sigma-Aldrich, St Louis, MO, USA). PCRs targeted the amplification of the small subunit ribosomal gene (18S rDNA), using both universal and myxozoan specific primers (Table 1). PCR reactions and cycling conditions were set-up according to Rocha et al. [40]. Amplification was confirmed by electrophoresis in a 1% agarose 1Â tris-acetate-EDTA buffer (TAE) gel stained with GreenSafe Premium (NZYTech, Lisbon, Portugal). PCR products were purified using the ExoFast method and sequenced directly using a BigDye Terminator v1.1 from the Applied Biosystems kit (Applied Biosystems, Carlsbad, CA, USA), and ABI 3700 DNA analyzer (Perkin-Elmer, Waltham, MA, USA; Applied Biosystems, Carlsbad, CA, USA; Stabvida, Oeiras, Portugal).

Sequence assembly, distance estimation, and phylogenetic analysis
The obtained forward and reverse sequence segments belonging to each sample were separately aligned using ClustalW in MEGA11 software [35,48]. For distance estimation, a dataset was generated that comprised all newly obtained sequences, plus all sequences available for other mugiliforminfecting Myxobolus, sphaeractinomyxon types, and the triactinomyxon type of Székely et al., 2007 (DQ473515).
This dataset included all sequences with highest similarity score to the newly obtained sequences (above the 80.0% threshold), as determined through BLASTn (Basic Local Alignment Search Tool). Alignments were performed using MAFFT version 7 online, and distance estimation was carried out in MEGA 11, with the p-distance model and all ambiguous positions removed for each sequence pair.
The dataset used for phylogenetic analyses also included all newly obtained sequences, plus 53 representatives of mugiliform-infecting Myxobolus and sphaeractinomyxon types, and the sequences of Myxobolus khaliji (KC711053), Myxobolus miyairii (KT001495) and Henneguya cynoscioni (JN017203) included as outgroup. Alignments of this dataset were also performed using MAFFT version 7 online. Bayesian inference analyses were conducted in MrBayes v.3.2.6 [44], with the general time reversible model with gamma-shaped rate variations across sites (Invgamma) (GTR + I + G). Posterior probabilities were calculated using the Markov chain Monte Carlo method, in which four chains ran simultaneously for 1 million generations, with burn-in set at 25%, and trees sampled every 500 generations. Maximum likelihood was performed in MEGA11, for 1,000 replicates, with the general time reversible substitution model with estimates of invariant sites and gamma distributed among site rate variation (GTR + I + G) chosen as the best suited model, based on the lowest score of the Bayesian Information Criterion (BIC) and corrected Akaike Information Criterion (AIC).

Myxozoan survey and overall prevalence of infection
The myxozoan survey conducted in this study revealed the presence of myxospores either dispersed in the tissue samples or contained within cysts. Of the 13 specimens of C. labrosus examined, 10 were infected (76.9%) with Myxobolus (Myxosporea, Myxobolidae); infections by other myxozoan genera were not observed. Most specimens were found to be parasitized by more than one Myxobolus species, with myxospores present in different organs (see Table 2). Highest species diversity and overall prevalence of infection occurred in the gills, followed by the peritoneum and fins. Co-infection was determined in the swim bladder of specimen#2 and in the fins of specimen#13, based on the analysis of ABI chromatograms, which included numerous overlapping peaks, and suspect base-calling errors in regions other than at the extremities. The species infecting the kidney of specimen#7 failed to be sequenced due to the low amount of parasite DNA ( Diagnosis: White spherical cysts, 1-2 mm in size, located between rays of caudal fin. Myxospores spherical to subspherical in valvular view and ellipsoidal in sutural view, with 6-10 sutural edge markings. Two large equally sized pyriform CTA CGG AAA CCT TGT TAC G ACT2f, ACT3f, ACT5f, ACT6f, ACT8f [50] polar capsules located side by side at anterior pole, occupying a little less than half myxospore body, each containing polar tubule coiled in 3-5 turns (Figs. 1C, 1D, 4B, Table 3). Site of infection: Caudal fin.
Site of infection: Gills.
Material deposited: Series of phototypes of the hapantotype, deposited together with a representative DNA sample in the Natural History and Science Museum of the University of Porto, Portugal, reference CIIMAR.2023.75.
Etymology: The specific epithet "chelonari" refers to the host genus.
Molecular data: One 18S rDNA sequence with 1,968 nucleotides and GenBank accession number OQ319218, representative of five identical sequences that were separately obtained from 5 cysts collected from a single infected specimen. Distance estimation revealed highest similarity to M. peritonaei n. sp. (98.1%) and M. invictus n. sp. (97.8%), with all other sequences included in the analysis displaying similarity values lower than 96.9%.  Table 3).
Site of infection: Gills.
Material deposited: Series of phototypes of the hapantotype, deposited together with a representative DNA sample in the Natural History and Science Museum of the University of Porto, Portugal, reference CIIMAR.2023.76.
Etymology: The specific epithet "aestuarium" is a noun referring to the estuarine habitat of the host species.
Molecular data: One 18S rDNA sequence with 1,972 nucleotides and GenBank accession number OQ319169, representative of two identical sequences that were separately obtained from both infected specimens. Distance estimation did not reveal similarity values higher than 96.6% to all other sequences reported here or retrieved from the NCBI database.
Site of infection: Gills.
Material deposited: Series of phototypes of the hapantotype, deposited together with a representative DNA sample in the Natural History and Science Museum of the University of Porto, Portugal, reference CIIMAR.2023.77.
Etymology: The specific epithet "invictus" refers to the city of Porto, historically known as "Invicta" (meaning: the invincible city).
Molecular data: One 18S rDNA sequence with 1,953 nucleotides and GenBank accession number OQ319224. Distance estimation did not reveal similarity values higher than 97.9% to all other sequences reported here or retrieved from the NCBI database.  Table 3).
Material deposited: Series of phototypes of the hapantotype, deposited together with a representative DNA sample in the Natural History and Science Museum of the University of Porto, Portugal, reference CIIMAR.2023.78.
Etymology: The specific epithet "douroensis" refers to the type locality, which is the Douro River.
Molecular data: One 18S rDNA sequence with 1,989 nucleotides and GenBank accession number OQ319165. Distance estimation did not reveal similarity values higher than 90.5% to all other sequences reported here or retrieved from the NCBI database.
Prevalence of infection: 1 infected of 13 specimens examined (7.7%) ( Etymology: The specific epithet "abdominalis" refers to the distribution of the parasite in the host abdominal cavity.
Molecular data: One 18S rDNA sequence with 1,982 nucleotides and GenBank accession number OQ319155, representative of two identical sequences that were separately obtained from both infected specimens. Distance estimation did not reveal similarity values higher than 93.3% to all other sequences reported here or retrieved from the NCBI database.
Site of infection: Gallbladder.
Material deposited: Series of phototypes of the hapantotype, deposited together with a representative DNA sample in the Natural History and Science Museum of the University of Porto, Portugal, reference CIIMAR.2023.82.
Etymology: The specific epithet "cucurbitiformis" refers to the pumpkin-shape of the myxospores.
Molecular data: One 18S rDNA sequence with 1,988 nucleotides and GenBank accession number OQ319220. Distance estimation matched the newly obtained sequence with that of Sphaeractinomyxon type 5 of Rangel et al. [29] (KU569314) (99.8%), with all other sequences included in the analysis displaying similarity values lower than 91.1%.
Material deposited: Series of phototypes of the hapantotype, deposited together with a representative DNA sample in the Natural History and Science Museum of the University of Porto, Portugal, reference CIIMAR.2023.83.
Etymology: The specific epithet "intestinicola" refers to the site of infection. Molecular data: One 18S rDNA sequence with 1,963 nucleotides and GenBank accession number OQ319222. Distance estimation matched the newly obtained sequence with that of Sphaeractinomyxon type 8 of Rangel et al. [29] (KU569317) (99.9%), with all other sequences included in the analysis displaying similarity values lower than 95.6%.

Differential diagnosis
Differentiation between the novel isolates reported in this study relied on the comparison of their 18S rDNA sequences, given that few significant morphometrically distinguishable features could be found between myxospores (see Table 3). Morphological comparisons were performed in relation to all Myxobolus spp. that were previously reported from mullet hosts and that are without molecular data.
All isolates reported in this study shared some morphometric similarity with Myxobolus hani Faye et al., 1999 (see Table 3), which could be readily differentiated based on different organ of infection (branchial spines) whenever compared to species infecting other organs, host genus (Mugil curema Valenciennes, 1836), geographic location (Atlantic Ocean off Senegal), and slight differences in terms of length and width [15]. Myxospores of M. hani have a slightly smaller width range than M. pinnula n. sp. and M. caudalis n. sp., and slightly greater length range than M. caudalis n. sp., M. invictus n. sp., M. douroensis n. sp. and M. labicola n. sp. In comparison to M. chelonari n. sp. and M. aestuarium n. sp., the myxospores of M. hani are thinner with a slightly greater length range, while being shorter than those of M. peritonaei n. sp., wider than those of M. abdominalis n. sp. and M. intestinicola n. sp., and having a greater width range than those of M. cucurbitiformis n. sp.
Except for M. aestuarium n. sp., M. abdominalis n. sp. and M. intestinicola n. sp., all other isolates further showed some morphometric similarity to Myxobolus cheni Schulman, 1962 (see Table 3), which could be promptly distinguished based on organ of infection (trunk muscle), host species [M. cephalus, and Planiliza haematocheilus (Temminck and Schlegel, 1845)], geographic location (China), and overall smaller myxospore width and longer polar capsules, also with slight differences in the length range for most cases [12].
Myxobolus peritonaei n. sp. further resembled Myxobolus parsi Das, 1996 (see Table 3) that, however, differs by having myxospores with a smaller length and width range, different organ of infection (gills), host species [Chelon parsia (Hamilton, 1822)], and geographic location (India) [9]. In turn, M. abdominalis n. sp. shared some morphometric similarity with Myxobolus adeli (Isjumova, 1964) Yurakhno and Ovcharenko, 2014 reported from several organs of Chelon auratus (Risso, 1810) in the Mediterranean Sea off Spain, Azov and Black Seas, and Myxobolus macrolepi Dorothy and Kalavati, 1992 from the intestine of Planiliza macrolepis (Smith, 1846) in India (see Table 3). Apart from these differences in organ of infection, host species and geographic location, the myxospores of M. adeli are larger with a greater thickness range and display three polar tubule coils instead of four, while those of M. macrolepsi are shorter with a greater width range, further displaying a higher number of polar tubule coils (6-7) [10,52].
Comparison with unnamed Myxobolus further revealed that most species described here share some morphometric similarity with the Myxobolus sp. 3 and Myxobolus sp. 5 reported by Rocha et al. [40] from C. ramada, with the first differing by having overall larger polar capsules with a higher number of polar tubule coils (6), and the latter differing by having unequally sized polar capsules. The myxospores of M. abdominalis n. sp. further resembled the Myxobolus sp. 1 and Myxobolus sp. 2 of Rocha et al. [40] with only slight morphometric variations, while those of M. intestinicola n. sp., which also shared some similar features with Myxobolus sp. 2 of Rocha et al. [40], could be differentiated based on the latter being less thick (see Table 3).
Overall, the comprehensive comparison of morphological, molecular, host-related, and geographic data performed in this study supported the description of these isolates as new species. The only exception was M. pupkoi, originally described from the gill arches of C. ramada in the Galilee Sea off Israel [21], and that was here identified from the gills of C. labrosus based on a 99.3% similarity score. This small genetic difference (0.7%) mostly took place at the 3 0 end of the 18S rDNA sequence, with a few nucleotide differences occurring elsewhere (position 694: T vs. A; positions 1508-1509: C vs. Y; and position 1538: A vs. W). Morphometric comparisons revealed significant differences among our isolates of M. pupkoi and the isolate reported by Gupta et al. [21] (see Table 3).

Phylogenetic analysis
Bayesian inference and maximum likelihood retrieved highly congruent tree topologies. Both analyses revealed the 18S rDNA sequences of the novel Myxobolus spp. clustering within Chelon-infecting lineages (Fig. 5). Myxobolus pinnula n. sp., M. caudalis n. sp., M. chelonari n. sp., M. aestuarium n. sp., M. invictus n. sp., M. peritonaei n. sp., M. abdominalis n. sp., M. intestinicola n. sp., and the novel sequence of M. pupkoi clustered within the subclades of the main Cheloninfecting lineage (Fig. 5, clade A1). The novel sequence of M. pupkoi specifically clustered with its conspecific sequence from the original species description (OL605966), while M. intestinicola n. sp. clustered together with its actinospore counterpart, the Sphaeractinomyxon type 8 of Rangel et al. [38] (KU569317). Myxobolus cucurbitiformis n. sp., M. douroensis n. sp., and M. labicola n. sp. clustered together with several sphaeractinomyxon types to form a novel independent, smaller Chelon-infecting lineage (Fig. 5, clade A2). Myxobolus cucurbitiformis n. sp. specifically clustered with the Sphaeractinomyxon type 5 of Rangel et al. [38] (KU569314), here determined as its actinospore counterpart. In accordance with previous cladograms [21,40], M. ramadus from the gills of C. ramada was retrieved at a more basal position (Fig. 5, clade A3), and is here shown as a sister branch to a small lineage of Myxobolus spp. that infect hosts of the genus Planiliza. A second lineage of Planiliza-infecting Myxobolus consistently branched at the basis of the most derived lineage of Chelon-infecting species. Finally, all Myxobolus spp. that infect hosts of the genus Mugil clustered together within the same lineage, while M. supamattayai formed a single branch as the sole representative of Crenimugil-infecting species (Fig. 5).

Discussion
Presently, it is widely accepted that the combination of several criteria, such as myxospore morphology and morphometry, host and tissue specificity/tropism, geographical location/ habitat, and molecular data, is required for the reliable description of myxozoan species [2,19]. Recent studies showed that, in the case of mugiliform-infecting Myxobolus, molecular data are mandatory for species identification, due to the notable lack of distinctive morphological characters that could allow myxospore differentiation [39,40]. Congruently, the Myxobolus spp. described in this study lacked significant morphologically distinguishable features allowing reliable differentiation amongst each other and in relation to other closely related species, thus reinforcing the need to perform sequencing of selected genetic markers whenever working with mugiliforminfecting Myxobolus. Morphology-based comparisons, however, remained necessary to establish differentiation from known mugiliform-infecting Myxobolus species that lack molecular data.
Worldwide, the thicklip grey mullet C. labrosus has been reported to host 5 Myxobolus spp. [39,40], including M. peritoneum and M. labrosus described from the peritoneum and urinary bladder of specimens caught from the Minho River estuary, in northern Portugal [40]. In this study, 11 new Myxobolus spp. are described from specimens obtained from the Douro River estuary, also located in northern Portugal. This species richness suggests that the real diversity of Myxobolus infecting C. labrosus is probably extremely underestimated. Although the Myxobolus community infecting this host varied between geographically close estuaries (the Douro and Minho River estuaries are less than 85 km apart), it should be considered that the study performed here was not ecologically significant given that seasonality was not evaluated. Therefore, a more thorough and organized effort in sampling C. labrosus specimens will most likely reveal even higher biodiversity of Myxobolus, with some species probably present in more than a single estuary. This hypothesis is supported by the matching of M. cucurbitiformis n. sp. and M. intestinicola n. sp. with Sphaeractinomyxon type 5 and Sphaeractinomyxon type 8, respectively reported by Rangel et al. [38] from the Aveiro estuary, 84 km south of the Douro estuary. In this study, the occurrence of highest species diversity and overall prevalence of infection in the gills agrees with previous studies that report this organ as a preferential site of infection for myxobolids, potentially allowing optimal myxospore release [14,31,32].
Even more patchy than our knowledge of the biodiversity of Myxobolus parasites, is that of their life cycles. The information presently available links this genus to a significant diversity of actinospore morphotypes [13]; however, a specific correspondence of mugiliform-infecting species to the sphaeractinomyxon collective group was recently proposed by Rocha et al. [42,43]. In this study, the molecular matching of M. cucurbitiformis n. sp. and M. intestinicola n. sp. with Sphaeractinomyxon type 5 and Sphaeractinomyxon type 8 of Rangel et al. [38], respectively reinforces the conjecture that sphaeractinomyxon types constitute specific life cycle counterparts of mugiliform-infecting Myxobolus. A 99.3% similarity score matched the sequence of one Myxobolus isolate from the gills of C. labrosus with that of M. pupkoi (OL605966) reported from the gill arches of C. ramada in the Galilee Sea off Israel [21], despite morphological comparisons showing significant differences between these geographical isolates (see Table 3). Although genetic similarity higher than 98% is usually considered to be representative of intraspecific variation among myxozoans [17], the boundary between myxobolid intraspecific and interspecific variability is not clear; there is no exact value of genetic similarity in the 18S rDNA sequences that determines whether two isolates can be discriminated as being conspecific or congener. Numerous instances of high values of intraspecific variation have been reported in myxobolids, namely M. koi (3.0%), H. corruscans  [7,11,16,24,33]. Since characters like morphology, host specificity, tissue tropism, and geographical location all come into play when inferring about intra-or interspecific variation, the decision of determining genetic differences as being representative of one or the other must be made individually. In the case of M. pupkoi, the small genetic difference determined between geographic isolates (0.7%) mostly took place at the 3 0 end of the 18S rDNA sequence and, therefore, cannot confidently be taken as a premise of interspecific variability. Furthermore, our isolate shares the same site of infection and host genus with the Israeli isolate, with Gupta et al. [21] even proposing that M. pupkoi was introduced into the Sea of Galilee via the Mediterranean Sea. Thus, the only considerable difference between our isolate and the Israeli isolate is myxospore morphometry, given that the myxospores observed in this study significantly surpassed the size range reported by Gupta et al. [21]. However, morphological plasticity is not uncommon between myxozoans, and has been reported between geographical isolates of several Myxobolus spp. [20,21,26,47,49,51]. For instance, in clear parallelism to our case, Israeli isolates of Myxobolus exiguus Thélohan, 1895 were also reported to have smaller myxospores than Portuguese isolates belonging to the same species, with identification based on the comparison of~2,000 bp sequences [21]. Thus, the genetic difference determined here between geographical isolates of M. pupkoi is proposed to be representative of intraspecific variability, with this species displaying morphological plasticity that may be hypothesized to correlate with adaptation to distinct environmental pressures and invertebrate communities.
In agreement with previous cladograms (e.g. [8,21,40]), our phylogenetic analysis showed all novel Myxobolus sequences clustering within the clade of mugiliform-infecting Myxobolus, which reinforces the proposed monophyletic origin of this group [40]. Congruently with host affinity being the main evolutionary driver of myxobolid radiation [6,27], hostassociated clustering was retrieved within the mugiliforminfecting clade, which included well-defined lineages of species parasitizing mullets belonging to the genera Chelon, Mugil, Planiliza, and Crenimugil. The identification of more than one Chelon-infecting lineage suggests that myxobolids acquired this fish genus as secondary hosts multiple times during their evolution. The same can be hypothesized for Planiliza-infecting species, which appear to be included within two separate lineages. Finally, the considerably high number of unmatched sphaeractinomyxon types included in the most derived Chelon-infecting lineages strengthens our contention that Myxobolus diversity in C. labrosus most likely remains underestimated. This shows a clear need to perform further studies targeting not only this species, but also others of the genus Chelon, and mullets in general.